#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <glpk.h>

using namespace std;

vector<double> Find_XI(vector<double> _beta, int i, double eps)
{
    vector<double> X;
    for(double j = _beta[i];j < _beta[i+1];j+=eps)
    {
        X.push_back(j);
    }
    return X;
}

int main()
{
	int index = 0, final_index = index;
	double final_ans = 0.0;
	for(double c = 0.8;c <= 0.9;c+=0.1) //0.710
	{
		double best_ans = 0.0, best_rho1=-1, best_rho2=-2;
		
		for(double rho1 = 0.65;rho1 <= 1.0;rho1 += 0.01) //0.245
		{
			for(double rho2 = 0.25;rho2 <= 1-rho1;rho2 += 0.01) // 0.718
			{

				ofstream file;
				file.open("LP.txt", ios::trunc);
				/**Param**/
				file << "param k;" << endl;
				file << "param rho1;" << endl;
				file << "param rho2;" << endl;
				file << "set I:= {0 .. k};" << endl;
				file << "set II:= {1 .. k};" << endl;
				file << "set III:= {0 .. k-1};" << endl;
				file << "param beta{I};" << endl;
				file << "param exp_beta{I};" << endl;
				file << "param exp_beta_rho1{I};" << endl; // exp(-beta[i] * rho1)
				file << "param exp_beta_rho2{I};" << endl; // exp(-beta[i] * rho2)
				file << "param exp_beta_oneminusrho{I};" << endl; // exp(-beta[i] * (1.0-rho1-rho2))
				file << "param exp_beta_rho{I};" << endl; // exp(-beta[i] * (rho1+rho2))
				int k = 2500;
				// double rho1 = 0.71000, rho2 = 0.24400, c = 0.71800;

				// double rho1 = 0.64863, rho2 = 1-rho1, c = 0.72941;
				vector<double> beta(k+1);
				for(int i=0;i<=k;i++)
				{
					beta[i] = c + 0.01 * i;
				}


				/**Var**/
				file << "var resultofthreephase;" << endl;
				file << "var Lambda1;" << endl;
				file << "var Lambda2;" << endl;
				file << "var Lambda3;" << endl;
				file << "var Delta{I};" << endl;
				file << "var delta{III};" << endl;
				file << "var v{I};" << endl;

				/**Objective Function**/
				file << "minimize z: resultofthreephase;" << endl;
				
				/**Constraint**/
				file << "subject to result3: resultofthreephase >= Lambda1 + Lambda2 + Lambda3;" << endl;
				file << "subject to Low_Lambda1: " << 1.0 - exp(-beta[0] * rho1) << " * (v[0] + Delta[0] / beta[0]) <= Lambda1;" << endl;
				file << "subject to Low_Lambda2: " << "sum {i in II} (exp_beta_rho1[i-1] - exp_beta_rho1[i]) * (1 - exp_beta_rho2[i]) * (beta[i-1] / beta[i] * v[i] + Delta[i-1]/beta[i]) <= Lambda2;" << endl;
				file << "subject to Low_Lambda3: " << "sum {i in II} (1 - exp_beta_oneminusrho[i]) * (beta[i-1] / beta[i] * v[i] + Delta[i-1]/beta[i]) * (" << exp(-rho1 * c) << "*(exp_beta_rho2[i-1] - exp_beta_rho2[i]) - rho2 / (rho1 + rho2) * (exp_beta_rho[i-1] - exp_beta_rho[i])) <= Lambda3;" << endl; 
				file << "subject to OPT: " << "1 <= v[k] + Delta[0] + sum {i in III} delta[i];" << endl;
				for(int i=0;i<k;i++)
				{
					file << "subject to Up_delta_one_C" << i << ": delta[" << i << "] <= " << 1.0 - exp(-beta[i+1]) << " * (v[" << i << "] - v[" << i+1 << "]);" << endl;
				}
				for(int i=0;i<k;i++)
				{
					file << "subject to Up_delta_two_C" << i << ": delta[" << i << "] <= Delta[" << i+1 << "] - Delta[" << i << "];" << endl; 
				}
				for(int i=0;i<k;i++)
				{
					vector<double> Xi = Find_XI(beta, i, 0.001);
					int index1 = 0;
					for(auto it=Xi.begin();it!=Xi.end();it++)
					{
						double gamma = (*it - beta[i]) / (exp(-*it) - exp(-beta[i+1]));
						file << "subject to Up_delta_three_C" << i << "_" << index1 << ": " << gamma << " * delta[" << i << "] <= Delta[" << i+1 << "] - Delta[" << i << "] - " << beta[i] - gamma * (1.0 - exp(-*it)) << " * (v[" << i << "] - v[" << i+1 << "]);" << endl;
						index1++;
					}
				}
				for(int i=0;i<k;i++)
				{
					file << "subject to V_C" << i << ": v[" << i+1 << "] <= v[" << i << "];" << endl;
				}
				file << "subject to V_C" << k << ": v[" << k << "] >= 0;" << endl;
				file << "subject to Low_Delta_C0: Delta[0] >= 0;" << endl;
				for(int i=0;i<k;i++)
				{
					file << "subject to Low_delta_C" << i << ": delta[" << i << "] >= 0;" << endl;
				}
				/*
				for(int i=0;i<=k-1;i++)
				{
					file << "subject to New_const_1" << i << ": delta[" << i << "] <= (v[" << i << "] - v[" << i+1 << "]) * " << 1.0 - exp(-beta[i+1]) << ";" << endl; 
				}
				
				for(int i=0;i<=k-1;i++)
				{
					file << "subject to New_const_2" << i << ": delta[" << i << "] <= 1.0 - " << 1.0 - exp(-c) << " * (v[0] + v[" << i+1 << "] - v[" << i << "]);" << endl;  
				}
				
				file << "subject to New_LowOPT: (v[k] * " << exp(-beta[k]) << ") + (v[0] * " << 1.0 - exp(-c) << ") + sum {i in III} v[i] * (exp_beta[i] - exp_beta[i+1]) >= 1 - Delta[0];" << endl;

				for(int i=0;i<=k-1;i++)
				{
					file << "subject to New_const_3_" << i << ": delta[" << i << "] <= " << exp(-beta[i]) << " * (Delta[" << i+1 << "] - Delta[" << i << "]);" << endl;
				}
				*/
				/***Data*/
				file << "data;" << endl;
				file << "param k := " << k << ";" << endl;
				file << "param rho1 := " << rho1 << ";" << endl;
				file << "param rho2 := " << rho2 << ";" << endl;
				
				file << "param beta :=" << endl;
				for(int i=0;i<=k;i++)
				{
					file << i << " " << beta[i] << endl;
				}
				file << ";" << endl;

				file << "param exp_beta_rho1 :=" << endl;
				for(int i=0;i<=k;i++)
				{
					file << i << " " << exp(-beta[i] * rho1) << endl;
				}
				file << ";" << endl;

				file << "param exp_beta_rho2 :=" << endl;
				for(int i=0;i<=k;i++)
				{
					file << i << " " << exp(-beta[i] * rho2) << endl;
				}
				file << ";" << endl;

				file << "param exp_beta_oneminusrho :=" << endl;
				for(int i=0;i<=k;i++)
				{
					file << i << " " << exp(-beta[i] * (1.0-rho1-rho2)) << endl;
				}
				file << ";" << endl;

				file << "param exp_beta_rho :=" << endl;
				for(int i=0;i<=k;i++)
				{
					file << i << " " << exp(-beta[i] * (rho1+rho2)) << endl;
				}
				file << ";" << endl;
				file << "param exp_beta :=" << endl;
				for(int i=0;i<=k;i++)
				{	
					file << i << " " << exp(-beta[i]) << endl;
				}
				file << ";" << endl;
				file.close();

				glp_prob *mip;
				glp_tran *tran;
				int ret;
				mip = glp_create_prob();
				tran = glp_mpl_alloc_wksp();
				ret = glp_mpl_read_model(tran, "LP.txt", 0);
				if(ret != 0)
				{
					fprintf(stderr, "Error on translating model\n");
					// goto skip;
				}
				ret = glp_mpl_generate(tran, NULL);
				if(ret != 0)
				{
					fprintf(stderr, "Error on generating model\n");
					// goto skip;
				}
				glp_mpl_build_prob(tran, mip);
				glp_scale_prob(mip, GLP_SF_AUTO);
				glp_smcp parm;
				glp_init_smcp(&parm);
				// glp_exact(mip, &parm);
				// parm.meth = GLP_DUALP;
				parm.presolve = GLP_ON;
				glp_simplex(mip, &parm);
				ret = glp_mpl_postsolve(tran, mip, GLP_SOL);
				if(ret != 0)
				{
					fprintf(stderr, "Error on postsolving model\n");
				}
				// glp_print_sol(mip, "123Result.txt");
				// ret = glp_print_sol(mip, NULL);
				double ans = glp_get_obj_val(mip);
				// final_ans = max(final_ans, ans);
				// index++;
				if(ans >= best_ans)
				{
					best_ans = ans;
					best_rho1 = rho1;
					best_rho2 = rho2;
				}
			//	if(final_ans == ans)
			//	{
			//		final_index = index;
			//	}
				glp_mpl_free_wksp(tran);
				glp_delete_prob(mip);
		/**	
			skip: glp_mpl_free_wksp(tran);
				  glp_delete_prob(mip);	
		**/
			}
		}
				ofstream file1;
				file1.open("Result_3Phase.txt", ios::app);
				file1 << "When c is " << c << ", the best rho1 is" << best_rho1 << ", rho2 is " << best_rho2 << ", ratio: " << best_ans << endl;
				file1.close();
	}
	// cout << "The Maximum ratio is " << final_ans << ", and the index is " << final_index << endl;
    return 0;
}
